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Abstract 

We investigate nonstationary excitations in 3D-Bose-Einstein condensates in a 
spherically symmetric trap potential under the modulation of scattering length 
with slowly varying frequencies (adiabatic modulation). By numerically solving the 
Gross-Pitaevskii equation we observe a step-wise increase in the amplitude of os- 
cillation due to successive phase locking between driving frequency and nonlinear 
frequency. Such a nonstationary excitation has been shown to exist by an analytic 
approach using variational procedure and perturbation theory in the action-angle 
variables. By using a canonical perturbation theory, we have identified the suc- 
cessive resonance excitations whenever the driven frequency matches the nonlinear 
frequency or its subharmonics. 
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1 Introduction 



The successful experimental realization of trapped Bose-Einstein condensates 
in alkali metal atoms has triggered immense interest in understanding the 
various properties of the ultra cold matter [Tf2] . The properties of the conden- 
sate wave function is usually described by a mean field Gross-Pitaevskii equa- 
tion [3] . For the past couple of years, there has been increased interest in study- 
ing the properties of Bose-Einstein condensates with time varying trap poten- 
tials and scattering lengths both experimentally and theoretically [H07IBIH] ■ 
In particular, temporal periodic modulation of scattering length by exploiting 
a Feshbach resonance is given a central importance in recent times. Earlier 
studies show that, in certain circumstances, periodically varying scattering 
length can stabilize the collapsing condensate piTfSfU] . However, very recently 
it has been shown that for a sign alternating nonlinearity an increase in the 
frequency of oscillations accelerates collapse [TO] . 

From another point of view, these periodic modulations lead to resonance 
phenomena in the condensates. Resonance is an interesting feature of an os- 
cillation under the action of an external periodic force manifesting in a large 
amplitude, when the frequency of the external force equals an integral mul- 
tiple of the natural frequency of oscillation. Although, the phenomenon of 
resonance is well understood in linear systems, the nonlinearity that arises in 
the Bose-Einstein condensates due to the inter-atomic interactions leads to an 
important problem of nonlinear resonance of wide interest. 

In general, there are two ways of driving a nonlinear oscillator by small driv- 
ing periodic force, namely, external and parametric. If the driving frequency 
is constant (that is, driving force is exactly periodic with time), the initial 
growth of the oscillator's amplitude with time is arrested by its nonlinearity. 
On the other hand, if the driving frequency is slowly varying with time (driv- 
ing force is almost periodic with time), the oscillator can stay phase locked 
but, on an average, increase its amplitude with time or a persistent growth in 
the amplitude takes place, and this phenomenon is known as autoresonance. 
In such systems, the nonlinear frequency is slowly varying with time. This 
autoresonance phenomena is also referred as adiabatic nonlinear phase lock- 
ing and synchronization. However, in certain nonlinear systems, frequency is 
independent of the amplitude of oscillation [TT1[T2"] . In these systems, when 
the driving frequency is slowly varying with time, one might expect succes- 
sive resonance excitations at subharmonic frequencies. We refer this type to 
resonance as a kind of subharmonic autoresonance [T5] . 

In the present work, we identify such a subharmonic autoresonance in Bose- 
Einstein condensates where it shows successive resonance excitations due to 
periodic modulation with slowly varying frequency. In Bose-Einstein conden- 
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sates it is experimentally feasible to vary the scattering length by either mag- 
netically or optically inducing a Feshbach resonance [UH]. Earlier works on 
Bose-Einstein condensation report the periodic modulation with constant fre- 
quency of scattering length [5117115115] . Along these lines, it is of potential in- 
terest to understand the dynamics of Bose-Einstein condensates under the 
action of periodically varying scattering length with slowly varying frequency. 
Motivated by the above, in this paper, we study the effect of such an exci- 
tation on the 3D Bose-Einstein condensates in a spherically symmetric trap 
potential and investigate the nature of parametric resonance. In particular, we 
point out the step-wise increase in amplitude of oscillation due to successive 
phase locking between the driving frequency and the nonlinear frequency (or 
its subharmonics) of the system. 

This paper is organized as follows. In Sec. [2] we discuss the properties of 
Bose-Einstein condensates driven by a periodic force with time varying fre- 
quency from numerical simulations using pseudo-spectral method. In Sec. [3] 
we describe the variational procedure and derive a reduced system of ordinary 
differential equations (ODEs) to describe the dynamics of condensate width. 
In sec. H] we analyze the width dynamics using a perturbed action angle vari- 
able theory for the reduced system of ODEs. Then in Sec. El by applying a 
canonical perturbation theory, we deduce the approximate nonlinear frequen- 
cies which are responsible for successive resonance excitation in BEC under 
the periodic modulation with slowly varying frequency. Finally in Sec. [61 we 
give a summary and conclusions. 



2 Nonlinear Gross-Pitaevskii Equation 



At ultra low temperatures, the time-dependent Bose-Einstein condensate wave 
function ^(r; r) at position r and time r may be described by the following 
mean-field nonlinear GP equation [3], 



^ + V{v) + gN\n^r)\ 2 -ih^- T 



0. 



Here m is the mass and N is the number of atoms in the condensate, g = 
4:%h 2 a/m is the strength of interatomic interaction, with a being the period- 
ically varying atomic scattering length. The normalization condition of the 
wave function is / dr|^(r; r)| 2 = 1. The three-dimensional trap potential is 
given by V(r) = |m \0J 2 x 2 + io 2 y 2 + oj 2 z 2 \ where u x , u> y , and u z are the an- 
gular frequencies in the x, y and z directions, respectively, and r = (x, y, z) is 
the radial vector (The standard notation (x,y,z) is reserved below for rescaled 
variables) . 
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2.1 Spherically symmetric trap potential 

In the spherically symmetric trap, i.e., u x = uj y = u z = ojq the trap potential is 
given by V(r) = \mio 2 f 2 , where u>q is the angular frequency and f the radial 
distance. The wave function can be written as ^(r;r) = tp{f,r). After a 
transformation of variables to dimensionless quantities defined by r = y/2f/l, 

t = tlu , I = ^(h/muo) and 0(r,t) = <p(r,t)/r = i/j(f, r)(4vr/ 3 ) 1/2 , the GP 
equation for the ground state wave function becomes 

<p(r,t) = 0, (2) 

where g = Na/l. The normalization condition for the wave function is then 

poo 

/ \<p(r,t)\ 2 dr = 2V2. (3) 
Jo 



2.2 Numerical results 

In order to study the nature of the condensate wave function, we solve the GP 
equation (j2J) numerically using a pseudo-spectral method [15]. In the pseudo- 
spectral method the condensate wave function is expanded in terms interpo- 
lating polynomials. When this expansion is substituted into the GP equation, 
the (space) differential operators operate on a set of known polynomials and 
generate a differentiation matrix operating on the unknown coefficients. Con- 
sequently, the time-dependent partial-differential nonlinear GP equation in 
space and time variables is reduced to a set of coupled ordinary differential 
equations in time which can then be solved by a fourth-order adaptive step-size 
controlled Runge-Kutta method [16] . 

A similar pseudo-spectral method has been used in [T7] with Hermite polyno- 
mials as the interpolant for the case of completely anisotropic trap potential. 
However, for the spherically symmetric condensates described by Eq. (j2J), it 
will be more advantageous to employ the integration in reduced dimensions 
with Laguerre polynomials as the interpolant. The Laguerre polynomials are 
more suited to the present problem as they are well defined in the interval 
r G [0, oo) and satisfy the boundary conditions of the wave function of the GP 
equation [T8] . 

We consider the case of periodically varying scattering length which can be 
achieved by means of Feshbach resonance. In order to take into account the 



d 2 r 2 
dr 2 4 ^ 



(p(r,t) 


2 d 


r 
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temporal modulation the scattering length can be taken as 



a(t) = e + ei cos 



uo(t) dt 



(4a) 



and 



g(t) 



Na(t) 



(4b) 



where 



u{t) = 2u + 5 — fit. 



(4c) 



is the frequency of the applied field with /i C 1 and 5 is a constant. We choose 
the parameter, without loss of generality and within the critical threshold limit 
for collapse, N — 1. The trap frequency ujq is fixed at 1. In order to study 
the resonance dynamics of (T5]), we calculate the root mean square distance, 

( r ) = \/v*)j which is defined as 



<r 2 > 



r 2 \ip(r, t)\ 2 dr. 



(5) 



When 5 = and fi = 0, the above system (TSJ) exhibits nonlinear resonance in 
which the steady growth of the (r) is arrested by the nonlinear frequency of 
the system |17j . Fig. [1] shows the plot of (r) as a function of time obtained 
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1000 1500 
t (in units of oJq 1 ) 



2500 



Fig. 1. Plot showing the rms value (r) as a function of time in the case constant 
frequency of the temporal modulation of the scattering length. The parameters are 
chosen as eo = 0, ei = 0.01, /i = and 5 = in Eq. (J4]). 

by solving (J2]) numerically for 5 = and /i = 0. Fig. [1] clearly indicates 
that the steady growth of (r) is controlled by the nonlinearity. However, for 
5 = and \i = 2 x 10 -3 , one observes an interesting resonance phenomenon 
giving rise to a step- wise increase in (r) as illustrated in Fig. [21 To understand 
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1000 2000 3000 4000 5000 

t (in units of uJq 1 ) 

Fig. 2. Plot showing the rms value (r) as a function of time in the case slowly varying 
frequency of the temporal modulation of the scattering length. The parameters are 
chosen as N = 1, loq = 1, e = 0, e\ = 0.01, fx = 2 x 10 -3 and 5 = in Eq. (|4j). 

such resonance excitations, we study the GP equation ([2]) by reducing it to 
a set of ODEs using a variational method and then by applying a canonical 
perturbation theory. 



3 Variational procedure 



The variational method is one of the simplest way to analyze the dynamics 
of Bose-Einstein condensate [7(191120] . In this method the original GP equa- 
tion (j2J) is reduced to a system of ODEs, with fewer variables describing the 
condensate wave packet, by assuming suitable trial wavefunction [7|19f20j . For 
convenience, Eq. (j2J) can be rewritten as (using the definition ip{r, t) = r0(r, t)) 



dt 



d 2 (f> 2d(f> | r 2 | 2 
dr 2 r dr 4 



(6) 



According to the variational method we assume the Gaussian wave function 
in the form ITll8ll9lTTOll211 



(r,t) = A(t) exp 



r2 



2a(t) 



(7) 



where A(t), a(t), b(t) and S(t) are amplitude, width, chirp and phase, respectively [20J, 
which are to be determined. 



The equation for the wave packet parameters A(t), a(t), b{t) and S(t) can be 
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obtained by calculating the average Lagrangian as [20 



where 



Lit) = 4tt / r 2 C(r,t)dr, 
Jo 



C(r,t) = - Ml + |V0| 2 + j|0| 2 + ^|0| 4 



(9) 



is the Lagrangian density. 



One can easily see that the above Gaussian form assumption (0) leads to the 
Lagrangian density, 



C{r,t) 



A *r 2 ( b - + - A +b 2 + -\+AH 
\ 2 a 4 4 / 



e _ r2/a2 + gify (10) 



and the average Lagrangian is given by 



L(t) 



3tt 3 / 2 



b 1 



^ 2 « 5 (S + 4 + & 2 + 7) + ^' 2 A 2 a y 5 + ^-g(t)A^. (11) 



3/2 



2 V 2 ' fl4 ' " ' A ) ' " ' 8^ 4 

The normalization condition for the wave function <p(r, t) as 

t)\ 2 dr = 2V2. (12) 



r 2 \0[r, 



From the above condition we get A 2 = Jj^- Therefore the average Lagrangian 
density can be rewritten as 

Lit) = 16 ^9(t) + 12v/2tt + 3v ^ 7rfl2 + 12 ^2 n b 2 a 2 + 6V2~7ra 2 b + 8V2n5. 
a 6 a 2 

(13) 

Then the Euler-Lagrangian equations for the Lagrangian L(t) lead to the 
following equations for the width a(t) and chirp b{t): 



a t = 2ab, b t 



Y^cr a 4 



(14) 



On eliminating b from the above Eq. (1141) . one is lead to the following evolution 
equation for a, 



4 8g{t) 
a tt + a = — + —=—-. 

a 6 A/7ra 4 

The above equation ffl~5]) can be rewritten as 



a tt + a = — 



4 , e(t) 



a 3 ' a 4 ' 



(15) 



(16) 
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where 



e(t) 



8N 



e + ei cos 



u(t) dt 



(17) 



uj(t) =2u + 5 — fit. 



For e(£) = 0, Eq. ffTB"]) is nothing but the well known Pinney equation |22j . 
This equation is completely solvable and the solution takes the following form 



a{t) = A + VA 2 -4sin(2t + 5) 



(19) 



Here A and 5 are arbitrary constants. The above nonlinear system Eq. (IT6|) 
with e(t) =0 possesses the interesting dynamical property that the frequency 
of the system is completely independent of the amplitude, unlike the case 
of standard nonlinear oscillators. For exceptions, see ref. [11]. However, for 
e(t) 7^ one has to study the full system Eq. f[T5"l) to understand the dynamics. 




t (in units of lu 



Fig. 3. Plot showing the width a as a function of time calculated from equation 
(Tl6|) . The parameters are N = 1, Uq = 1, eo = 0, e% = 0.01, (jl = 2 x 10 -3 and 
(5 = 0. 



In order to study the nature of oscillation for the width a we solve the above 
equation f)16p numerically using a fourth order Runge-Kutta method with 
appropriate initial conditions. For /x ^ (/i = 0.002), we observe that a step- 
wise increase in the amplitude of oscillation in a as shown in Fig. [3j However, 
when jj, = this shows a usual nonlinear resonance phenomena [20J . 

One can easily relate the width a(t) of the wavefunction obtaining from the the 
solution of the ODE (16) to the root mean square distance of the condensate 
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(r) as 



2^/2 V 



r 4 |(/)(r,i)| 2 (ir 



which can be evaluated using the from of the wavefunction <f)(r,t) [eq. (7)] 
and the fact that the amplitude of the wavefunction is given by the relation 
8 \/2 

A 2 = — — — obtained from the normalization condition (12). Thus the width 



of the condensate can be related to the root mean square distance of the 
condensate as (r) = a a 



4 Perturbation theory in the action-angle variables 



In order to have a more clear understanding on the resonance phenomena, 
we analyze the Eq. (fTEj) by constructing the action- angle variables. Treating 
Eq. (fT6|) as nearly integrable system(e(t) <C 1), it is more convenient to con- 
sider the action-angle variables in which one can employ a perturbation theory 



4-1 The action-angle variables for the unperturbed system 



First, let us consider the unperturbed system e(t) = 0. The Hamiltonian for 
the unperturbed problem is given by 

H = £ + U{a), (20) 



where 



TT, \ a2 2 



If E is the energy then the width, a(t), oscillates between a m i n = y E — \fE 2 — 4 



and a max = y E + E 2 — 4. Since the Hamiltonian is conserved, i.e., Hq = E, 

da 



V2E -a 2 - 4a- 2 
On integrating the above equation, we obtain 



dt. (21) 



(t) = y l E + VE 2 -4 sin 6, (22) 
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where 9 = 9q + 2t and 9q being integration constant, see Eq. (|T9l) . 



We now make a change to action-angle variables. In this problem, where the 
phase space is two-dimensional, there is one action variable / and a conjugate 
angle variable 9. The action variable is given by 

/ = — $> dfda 
2tt J 




2E - a 2 - 4 da. (23) 
a 2 

This integral can be easily evaluated by using the energy (E) to express the 
momentum at in terms of position a. The resulting integral is elementary and 
leads to 

I=\{E-2). (24) 

By expressing the total energy E in terms of the action J, the unperturbed 
Hamiltonian is written as 

H = E = 2(1 + 1), (25) 

which is a function of / only. The change of variables to action-angle variable 
is canonical, so that the Hamiltonian's equations retain their form 



d9 dH (I) 



dt dl 

where Qq = 2. 



n Q , (26b) 



4-2 The action-angle variables for the system with perturbation 



Now we consider the problem with perturbation. Then the Hamiltonian for 
perturbed problem fflBI) can be written as 

H = H + e(t)Ht = H + e(t)-L (27) 



where 



= { — — i cos(2t + 5t-/it 2 ), and e = 0. 
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Further, one can express H in terms of action angle variables I and 6 using 
the Eqs. (122]) and (125} . 



H(I,6) = H (I) + e(t)H 1 (I,ff), 



where 



i 

372' 



6V2 (l + / + y/I{2 + I) sin 
Thus the Hamiltonian equations of motion are given by 



(2f 



(29) 




1000 2000 3000 

t (in units of uJq 1 ) 



5000 



1.5 




1000 2000 3000 

t (in units of uJq 1 ) 



5000 



Fig. 4. Plot showing (a) the action / as a function of time obtained by solving the 
reduced system Eq. fl30f) and (b) the width a(t) obtained from eq. [31] as a function 
of t. The parameters are N = 1, ojq = 1, eo = 0, e\ = 0.01, /i = 2x 10 -3 and (5 = 0. 



Tt = ~ e{t) ^ 



(30a) 
(30b) 
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where 



dHi 

~~de~ 



1(2 + 1) cos 9 



sin 9 



5/2 ' 



dHi 
~~df 



1 + 



1 + / 

7(2 + 7) 



sin 9 



Ay/2 ( 1 + I + J7(2 + 7) sin# 



5/2 " 



The values of 7 and # are obtained, as a function of t, by solving Eqs. f[3"Ul) 
numerically using a fourth order Runge-Kutta method with appropriate initial 
conditions. Fig. Hl(a) shows the plot of 7 as a function of time. It is evident that 
I retains relatively small values. Substituting I and 9 values in equations (I25I) 
and (1221) . we get the expression for condensate width a in terms of action-angle 
variables as 



2(1 +7) +2JI(2 + 7) sin9 



(31) 



We note that the variation of a [in Eq. (l3T]) ] depicted in Fig. H](b), after using 
the results for I, compares well with that shown in Fig. [3j Further we have also 
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1000 2000 3000 4000 5000 



t (in units of uj ) 

Fig. 5. Plot of the (r) 

max obtained by solving eq. ([2]), from the numerical solutions 
of eq. (fT5]) and from perturbation theory ([30l 

compared, in Fig. [51 the results obtained from the direct numerical solution 
of the GP equation ([2]), using the numerical solution of (fl5|) obtained from 
variational procedure and the results of the perturbation theory by solving 
equations ( 1301) . 
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5 Canonical perturbation Theory 



In order to understand the resonance phenomena one has to find out the 
approximate nonlinear frequency of the above system which is responsible for 
the resonant excitations. For this purpose we adopt a canonical perturbation 
theory to find the nonlinear frequency of near integrable systems [23l24|l25j . 
The basic idea of this method is to rewrite the perturbed system (1281) into a 
new set of action angle variables (J, <p) by a canonical transformation to a new 
Hamiltonian K(J), which depends on the new action J only, that is, 

H(I,9)^K(J). 

For / C 1, one can express the perturbed Hamiltonian Hi in the binomial 
expansion form as, 



Hi{i,e) 



6V2 



1- - (7 + ^/(2 + /) sin# 



(32) 



We now attempt a canonical transformation to action-angle variables (J,</>) 
for H(I, 6) through a type II generating function [56] as 

3 (J, 6) = S (J, 9) + e(*)5i(J, 6) + e(t) 2 S 2 (J, 6) + ..., (33) 

where Sq = JO is the identity generator. Then the new action-angle variables 
are given by [26] 

' m (34a) 



89 
dS(J : 



dJ 



(34b) 



Rewriting the Hamiltonian (1281) in terms of J using the relation (134|) . such a 
transformation gives the new Hamiltonian K(J), that is, 

Ha (pmy e(t)Hl {^m, e y K(J) . (35) 

Furthermore, one now expands the new Hamiltonian K(J) as a power series 
of e(t), that is, 

K{J) = K {J) + e{t)K x {J) + t{t) 2 K 2 {J) + ... (36) 

By making a Taylor's series expansion of H and Hi in Eq. (|35|) . using the 
relation for S(J,9) from ([3*31 . and equating the coefficients of various powers 
of e(t), one can obtain a system of equations. The new Hamiltonian K(J) can 
then be obtained by solving this set of equations. 
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The zeroth order term Kq(J) = Hq of the new Hamiltonian is obtained by 
equating the coefficient of e(t)°. It is easy to see that K (J) can be simply 
found by replacing / by J in the zeroth order Hamiltonian H . Therefore the 
zero-order frequency is given by, 

^0 = ^ = 2 (37) 
Similarly, equating the coefficient of e(t), one obtains 

= Q Q -j£-H 1 (J,9) (38) 

At this point, one can exploit the periodicity of the motion in the angle variable 
9. Since Si, i = 1,2, . . ., are assumed to be periodic in 9. The averaging of 
Eq. (1381) over 9 leads to the mean of the derivatives of Si vanish. Then the 
first-order Hamiltonian correction is given by 

Kl (J)=l-~J. (39) 
Hence the first order correction to frequency is given by 



^ (J) 



dJ 
3 



(40) 



from Eq. (134"]) upto first order, 



J = I - <t) mJi, (41a) 

We get / and 9 values by solving Eq. (I3"U|) numerically and substituting these 
values in Eq. fj4T|) . we get the J value. 



The total nonlinear system frequency is, 

n(J) = tt (J) + e(t)^i(J). 

In the above expression the second term in the right hand side can be neglected 
as e(t) is very small and one can find Q(J) ~ ^o(^)- Here the nonlinear 
frequency is almost constant while the driving frequency is slowly varying 
with time. Hence one could expect a phase locking when the driving frequency 
matches with the nonlinear frequency or its integer multiples. 
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Fig. 6. Phase locking between nonlinear system frequency and driven frequency, 
(top) Plot of the (r) max obtained by solving Eq. ([2]) and (bottom) the driving fre- 
quency and the different harmonics of nonlinear frequency is shown. 

From Fig. [6] it is easy to see that at t ~ 1000, the driving frequency ui(t) 
matches with the nonlinear system frequency O (i.e, u(t) = — O), at which 
the first resonance jump occurs in the original problem as shown in Fig. 
It explains the primary resonance that occurs at t ~ 1000 in the original 
problem Eq.(j2J). At another time, t w 1500, where u(t) = —20, a subharmonic 
resonance occurs leading to a second jump in the amplitude as shown in Fig. [2j 
As a result a step-wise increase in the amplitude as shown in Fig. [2] due to 
the phase locking between drive and system whenever the driving frequency 
matches with the system nonlinear frequency or its integral multiples (ie, 
to(t) = nO, where n = 0, ±1, ±2, . . .). We refer this phenomena as a kind of 
subharmonic autoresonance. 



6 Summary and conclusions 



In summary, we have studied the autoresonant excitations in Bose-Einstein 
condensates under the action of external periodic modulation with time de- 
pendent frequency. By numerically solving the corresponding Gross-Pit aevskii 
equation we have observed that there occurs a successive phase locking with 
step-wise increase in the overall amplitude of oscillation. We have employed 
a variational procedure using Gaussian trial wave function to simplify the 
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problem in fewer coordinates. Then the reduced system has been found to 
show a kind of subharmonic autoresonance phenomenon with successive res- 
onance excitations. Further, we have studied the problem in the action-angle 
variables which allows one to make a perturbation analysis and thereby iden- 
tify the approximate frequency that is responsible for the successive resonant 
(step-wise) excitations. The results obtained from the canonical perturbation 
theory compare well with the numerical solution of the original problem. 
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